Dynamics in the presence of attractive patchy interactions 
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We report extensive monte-carlo and event-driven molecular dynamics simulations of a liquid 
composed by particles interacting via hard-sphere interactions complemented by four tetrahedrally 
coordinated short-range attractive ("sticky") spots, a model introduced several years ago by Kolafa 
and Nezbeda [J. Kolafa and L Nezbeda, Mol. Phys. 16187(1987)]. To access the dynamic properties 
of the model we introduce and implement a new event-driven molecular dynamics algorithm suited 
to study the evolution of hard bodies interacting, beside the repulsive hard-core, with a short- 
ranged inter-patch square well potential. We evaluate the thermodynamic properties of the model 
in deep supercooled states, where the bond network is fully developed, providing evidence of density 
anomalies. We show that, differing from models of spherically symmetric interacting particles, in 
a wide region of packing fractions (j) the liquid can be super-cooled without encountering the gas- 
liquid spinodal. In particular, we suggest that there is one optimal <j) (not very different from the 
hexagonal ice (j)) at which the bond tetrahedral network fully develops. We find evidence of the 
dynamic anomalies characterizing network forming liquids. Indeed, around the optimal network 
packing, dynamics fasten both on increasing and decreasing (j). Finally we locate the shape of the 
isodiffusivity lines in the {(j) ~ T) plane and establish the shape of the dynamic arrest line in the 
phase diagram of the model. Results are discussed in connection to colloidal dispersions of sticky 
particles and gel forming proteins and their ability to form dynamically arrested states. 
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I. INTRODUCTION 

This article presents a detailed numerical study of the 
thermodynamics and of the dynamics of a model intro- 
duced several years ago by Kolafa and Nezbeda Q as a 
primitive model for water (PMW). The model envisions 
a water molecule as a hard sphere (HS) whose surface is 
decorated by four short ranged "sticky" spots, arranged 
according to a tetrahedral geometry, two of which mimic 
the protons and two the lone-pairs. Despite its original 
motivation, the Kolafa and Nezbeda model is represen- 
tative of the larger class of particles interacting via lo- 
calized and directional interactions, a class of systems 
which includes, besides network forming molecular sys- 
tems, also proteins 2, 3, 4] and newly designed colloidal 
particles^. Indeed, recent developments in colloidal sci- 
ence are starting to provide particles with specific di- 
rectional interactions|a|. In the same way as sterically 
stabilized colloids have become the ideal experimental 
model for realizing the hard-sphere fluid, novel physical 
chemical techniques will soon make available to the com- 
munity colloidal analogs of several molecular systems. A 
colloidal water is probably not far from being realized. 

Recent work 7J has focused on the dynamics of col- 
loidal particles interacting with a restricted number of 
nearest neighbors. In Ref. (31, |3] particles are interacting 
via a limited- valency square well model0, 0, 0|, im- 
posing a many body constraint on the maximum num- 
ber Umax of bonded interactions. It has been found that 
when Umax < 6, a significant shrinking of the liquid-gas 
(or colloidal rich-colloidal poor) spinodal takes place. A 
window of packing fractions values opens up in which it is 
possible to reach very low temperature (and hence states 
with extremely long bond lifetime) without encounter- 



ing phase separation. This favors the establishment of 
a spanning network of long-living bonds, which in the 
colloidal community provides indication of gel formation 
but which, in the field of network forming liquids, would 
be rather classified as glass formation. The study of the 
dynamics of the PMW provides a test of the rimax = 4 
results, in the absence of many-body interactions and in 
the presence of a geometric correlation between the bond- 
ing sites, retaining the maximum valency. This article, 
by reporting results on a model which can be at the same 
time considered as a simple model for the new generation 
of patchy colloids or for network forming liquids, starts 
to bridge the gap between these two fields. 

Thermodynamic and structural properties of several 
primitive models for water (and other bonded sys- 
tems) have been studied in detail during the last 30 
years 0,0, 0,0, 01, since this type of primitive models 
have become one of the landmarks for testing theories of 
association^^ 0, 

theory of Wertheimjl^ ll7| has been carefully compared 
to early numerical studies, suggesting a good agreement 
between theoretical predictions and numerical data, in 
the temperature and packing fraction regions where it 
was possible to achieve numerical equilibration IS]. Re- 
cently, the increased numerical facilities, extending the 
range of studied state points, have clarified that de- 
viations from the theoretical predictions start to take 
place as soon as the number of bonds (between different 
patches) per molecule increases and a network of bonded 
particles appears poll2^ . Geometric correlations between 
different bonds, not included in the theory, are responsi- 
ble for the break down of the theoretical and numerical 
agreement. Attempts to extend the perturbation theory 
beyond first order do not appear to be able to cure the 
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problem[23|- The PMW is a good candidate for test- 
ing new theories of association and, for this reason, it 
is important to clearly establish numerically the low T 
behavior of the supercooled liquid state. The equilib- 
rium PMW phase diagram, recently calculated 15j, in- 
cludes two crystal regions and a metastable fluid-gas co- 
existence. 

All previous studies of primitive models for sticky di- 
rectional interactions have focused on thermodynamic 
and static properties of the model. But the ability of fully 
exploiting the fast developments taking place in colloidal 
physics i24, 25] requires understanding not only the equi- 
librium phases of systems of patchy particles and their 
modifications with the external fields, but also under- 
standing the kinetic phase diagram|2^, i.e. the regions 
in phase space where disordered arrested states can be 
expected, and when and how these states are kinetically 
stabilized with respect to the ordered lowest free energy 
phases. In this respect, it is worth starting to estab- 
lish the dynamic properties of simple models of patchy 
interactions, since the simplicity of these models (based 
on hard sphere and square well interactions) have the 
potentiality to provide us with an important reference 
frame and may play a relevant role in deepening our 
understanding of the dynamic arrest in network form- 
ing liquids, in connecting arrest phenomena associated 
to gel formation|3, 113 (the establishment of a percolat- 
ing network of long lived bonds) and arrest related to 
excluded volume effects and the dependence of the gen- 
eral dynamic and thermodynamic features on the number 
and spatial location of patchy interactions. The case of 
the PMW reported here is a good starting one. In this 
article we report thermodynamic data, extending the pre- 
viously available information to lower temperatures and, 
for the first time, dynamic information obtained solving 
the Newton equations using a new algorithm based on 
event-driven propagation. 

II. THE MODEL AND NUMERICAL DETAILS 

In the PMW, each particles is composed of an hard 
sphere of diameter a (defining the length scale) and by 
four additional sites located along the direction of a tetra- 
hedral geometry. Two of the sites (the proton sites H) 
are located on the surface of the hard sphere, i.e. at dis- 
tance 0.5(7 from the center. The two remaining sites (the 
lonc-pair sites LP) are located at distance 0.45cr. Besides 
the hard-sphere interaction, preventing different particles 
to sample distances smaller than a, only the H and LP 
sites of distinct particles interact via a square well (SW) 
potential usw of width S — 0.15a and depth mq, i.e. 

usw = -wq r <d (1) 
= r>S, 

where r is here the distance between H and LP sites. 
The choice of 6 = 0.15a guarantees that multiple bond- 
ing can not take place at the same site. The depth of the 



square well potential uq defines the energy scale. Bond- 
ing between different particles is thus possible only for 
specific orientations and distances. In the linear geom- 
etry, the maximum center-to-center distance at which 
bonding is possible is l.lcr since the LP site is buried 
0.05f7 within the hard-core, a value typical of short-range 
colloid-colloid interactions. 

We have studied a system of = 350 particles with 
periodic boundary conditions in a wide range of packing 
fraction (f) = ir/Gna^ (where n is the number density) 
and temperatures T, where T is measured in units of 
uq {ks = !)• We perform both Monte Carlo (MC) and 
event driven molecular dynamics. In one MC step, an 
attempt to move each particle is performed. A move 
is defined as a displacement in each direction of a ran- 
dom quantity distributed uniformly between ±0.05 a and 
a rotation around a random axis of random angle dis- 
tributed uniformly between ±0.5 radiant. Equilibration 
was performed with MC, and monitored via the evolution 
of the potential energy (a direct measure of the number 
of bonds in the system) . The mean square displacement 
(MSD) was also calculated to guarantee that each par- 
ticle has diffused in average more than its diameter. In 
evaluating the MSD we have taken care of subtracting the 
center of mass displacement, an important correction in 
the low T long MC calculations. At low T simulations re- 
quired more than 10^ MC steps, corresponding to several 
months of CPU time. 

We have also performed event driven (ED) molecular 
dynamic simulations of the same system, modeling par- 
ticles as constant density spheres of diameter a and mass 
TO. The momentum of inertia is diagonal and equal to 
ma^/ 10 . The algorithm implemented to propagate the 
newtonian trajectory in the presence of patchy square 
well interaction is described in details in Appendix I VIII 
In ED dynamics, time is measured in units of cta/to/uq- 
Assuming as to the mass of the water molecule, as uq 
a typical value for hydrogen bond (~ 20kJ/mol) and as 
a the nearest neighbor distance in water (0.28nTO), the 
unit of time corresponds ~ 0.3ps. All static quantities 
have been evaluated with both MC and MD configura- 
tions finding no differences. 

Pressure, measured in units of ug/a^, has been cal- 
culate as sum of three contributions. A trivial kinetic 
contribution, equal to nksT. A positive HS contribution 
and a negative contribution arising from the SW interac- 
tion. Details of the calculation of P in both MC and ED 
simulations is provided in the Appendix I Villi 



III. RESULTS: STATIC 
A. Potential Energy E 

Since in the PMW each site can take part to only one 
bond, due to geometric constraints fixed by the small 
value of i5, the lowest energy configuration is defined by 
four bonds per particles, corresponding to a ground state 
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energy per particle Eg 



-2 (in units of uq). Of course, 



this absolute ground state value may not be accessible 
at all (j), due to the strong constraints introduced by the 
bonding geometry. According to Wertheim's first order 
thermodynamic perturbation theory, the T and </> de- 
pendence of the potential energy per particle E is given 
byEmilll 



where 



0.5 



J 



E — Egs — 



1 + 192(e- - l)(f>J 



0.5 



Ci(l-0/2)-C20(l + ^ 
(1 ~ 0)3 



(2) 



(3) 



(4) 



with ci = 2.375 X 10-^ and ca = 2.820-**[T3. The 
Wertheim theory, which assumes uncorrelated indepen- 
dent bonds, predicts as low T limit of Eq.[21an Arrhcnius 
T-dependence, 



lim E-Egs = 



,-0.5/T 



V1920J 



(5) 



i.e. with an activation energy of half bond energy. It 
is worth observing that such an Arrhcnius law, with an 
activation energy equal to 0.5 characterizes the low T de- 
pendence of the energy in the n„iax modelf^ 1^ [a model 
of particles interacting via a SW potential with an ad- 
ditional constraint on the maximum number of bonds], 
where no geometric correlation between bonds is im- 
posed. 

Fig- n shows the T dependence of the potential energy 
for different isochores. As discussed in details in the fol- 
lowing, for (j) < 0.24 a phase separation is encountered on 
cooling, preventing the possibility of equilibrating one- 
phase states below T « 0.11. For cf) > 0.24 the system 
remains homogeneous down to the lowest investigated 
T. The low T-behavior is expanded in Fig. ^(bottom). 
With the extremely long equilibration runs performed, 
proper equilibration is reached only for T > 0.09. The 
enlargement of the low T region shows that the absolute 
ground state value — 2uo is closely approached at </> « 0.3. 
At higher or smaller 0, the potential energy appear to 
approach a constant value larger than —2uq. Consis- 
tent with previous claims^^, high T data are very well 
represented by first order thermodynamic perturbation 
theory. Systematic deviations between theory and simu- 
lation data appears as soon as the number of bonds per 
particle becomes bigger than one. Comparing the sim- 
ulation data with the Wertheim theory, it is confirmed 
that the physics of the network formation is completely 
missing in the first-order perturbation theory. 

Fig. 121 shows the p dependence of E along isochores. 
At high T {T > 0.13), a monotonic decrease of E is 
observed, caused by the increased bonding probability 




FIG. 1: Potential energy for the PMW. The top panel shows 
data for all studied isochores as a function of T. The lower 
panel shows an enlargement of the low T region, where the 
network is fully developed. Note that for this model, the 



lowest possible energy is Eg 



-2. 



induced by packing. In this T region, the number of 
bonds is at most of the order of two per particle. Com- 
pletely different is the situation for lower T. The (j) de- 
pendence becomes non-monotonic. There is a specific 
value of the packing fraction (0 0.3) at which the low- 
est energy states are sampled. In the following we define 
the optimal network packing fractions as the range of 
packing fractions for which it is possible to fully satisfy 
the bonds in a disordered homogeneous structure. At 
(f) « 0.3, the number of bonds at the lowest investigated 
T (the lowest T at which equilibration is feasible with 
several months of computation time) is about 3.8 per 
particle, i.e. about 95% of the bonds are satisfied. The 
range of optimal 0s appears to be rather small. Indeed 
for packing fractions lower or higher than this optimal 
(f) K, 0.314, the formation of a fully connected network 
is hampered by geometric constraints: at lower 0, the 
large inter-particle distance acts against the possibility 
of forming a fully connected network, while at large </>, 
packing constraints, promoting close packing configura- 
tions are inconsistent with the tetrahedral bonding ge- 
ometry. Not surprisingly, = 0.314 is within the range 
of (j) values which allow for a stable open diamond crystal 
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FIG. 2: Potential energy vs. (p along isotherms. Symbols: 
simulation data. Lines: Wertheim's theory. 



phase (0.255 < cf) < 0.34) [ly. A reduction of the geomet- 
ric constraints (as in the Umax modelQiH^) increases the 
range of optimal (p. It is worth also noting that the liquid 
side of the spinodal curve is close to the region of optimal 
network <f>. 

The existence of a convex form for the potential energy 
(here for (p ^ 0.3) has been observed in several other mod- 
els for tetrahedral networks, including models for water 
(and water itself|2^). It has been pointed out that a neg- 
atively convex (j) dependence is indicative of a destabiliza- 
tion of the free energy jl^ and a precursor of a possible 
liquid-liquid critical point (in addition to the lower (j) gas- 
liquid one). Liquid-liquid critical points have been ob- 
served in several models for water [3(1 l3ll l3^ l3^ 133. l35l| . 
Indeed, the Helmholtz free energy A is related to U (the 
sum of the kinetic and potential energy) via A = U — TS, 
where S is the entropy. The curvature of an isotherm of 
A must be positive for a homogeneous phase of a specified 
volume V to be thermodynamically stable. The curva- 
ture of A can be expressed as 



\dvy. 



-T 



(6) 



Since P = — (f^)^ the inverse compressibility Kt 
— l/V{dV/dP)T is related to the curvature of A by 



1 



= V 



dHj_ 



T 



(7) 



The curvature of A is thus proportional to 1/ Kt for fixed 
V. Since 1/ Kt must be positive for a thermodynamically 

stable state, for the range of V in which (^f^j!-^ < 0, the 

contribution of the internal energy reduces the thermo- 
dynamic stability of the liquid phase. The liquid remains 
stable where U has negative curvature only because the 
contribution of the entropic term in Eq. (Hjis large enough 
to dominate. Yet entropic contributions to these thermo- 
dynamic quantities are suppressed as T decreases, due to 
the occurrence of the factor of T in the second term on 



the right-hand side of Eq. |H1 Hence the U — V data sug- 
gest that at lower T a single homogeneous phase of the 
liquid will not be stable for certain values of V, leading to 
a separation into two distinct liquid phases of higher and 
lower volume. Due to the predominant role of E in the 
free energy at low T, the possibility of a phase separation 
of the PMW liquid into two liquid phases of different 0, 
for (j) > 0.3 and T lower than the one we are currently 
able to equilibrate should be considered. 




FIG. 3: Arrhenius representation {E — Egs vs 1/T) of the 
potential energy around the optimal network density. The 
dashed line, shown as a reference, has an activation energy of 

Fig.Olshows lii{E-Egs) vs 1/T. At the optimal (p, the 
energy of the fully connected state is approached with an 
Arrhenius law, characterized by an activation energy of « 
Iuq, clearly different from the 0.5 value predicted by the 
Wertheim theory. For larger (j) values, data suggest that 
the lowest reachable state has an energy different from 
—2uo, consistent with the expectation that on increasing 
(j), geometric constraints forbid the development of a fully 
connected network even at the lowest T. 



B. 



The Wertheim's prediction for the T and (j) dependence 
of the PMW pressure (the equation of state) is 

P = Phs (8) 
96(ei/^ - 1) ci0(l + (/) - 0.502) _ 2c2(^2(^ _^ 20) 



-nkgT 



{1 + cy 



(1 



where Phs is the pressure of the HS fluid at the same 
packing fraction. Phs is very well represented by the 
Carnahan-Starling E0S|3| 



Phs = nksT 



(1 + + 0^-03) 

(1 - 0)3 



(9) 



The Wertheim EOS predicts a vapor-liquid critical point 
at Tc = 0.1031 and 0c = 0.085ll5j. The vapor-liquid 
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spinodals calculated according to the Wertheim theory 
and from simulation data are reported in Fig. ^ The 
numerical estimate is provided by locating, along iso- 
chores, the highest T state point in which phase sepa- 
ration is observed and the T at which the small q limit 
of the structure factor is smaller than five. These two 
state points bracket the spinodal locus. It is interesting 
to compare the liquid-gas spinodal of the PMW with the 
corresponding spinodal of the symmetric spherical square 
well potential with same depth and well width i5 = 0.15. 
In that case, the critical point is located at « 0.56 and 
(pc ~ 0. 212(33 ^'^d high packing fraction (the liquid) 
side of the spinodal extends beyond — 0.6. The net re- 
sult of decreasing the surface available to bonding and of 
limiting to four the maximum number of nearest neigh- 
bors which can form bonds is the opening of a wide region 
of 4) values where (in the absence of crystallization) an 
homogeneous fluid phase is stable (or metastable). This 
finding is in full agreement with the recent work of 0, 
where a saturated square well model was studied for dif- 
ferent values of the maximum valency. Indeed, it was 
found that when the number of bonds becomes less then 
six, the unstable region (the surface in the (0 — T) plane 
encompassed by the spinodal line) significantly shrinks, 
making it possible to access low T states under single 
phase conditions. 
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FIG. 4: Thermodynamic phase diagram for the PMW. 
The theoretical Wertheim prediction for the locus at which 
{dP/dV)T = is compared with numerical estimates of the 
spinodal, calculated by bracketing it via the locus of point at 
which S{0) « 5 and the locus of points where a clear phase 
separation is detected. The location of the bond percolation 
line is also reported. 

Fig. [S] shows P{T) for different isochores. In agree- 
ment with previous analysis, P is well represented by the 
Wertheim theory only at high temperature. At low T sev- 
eral interesting features are observed: (i) for < 0.25, 
isochores end in the spinodal line, (ii) in the simula- 
tion data, a clear difference in the low T behavior is 




FIG. 5: Isochores of P according to the Wertheim theory 
(top) and as calculated from the simulation data (bottom). 
Symbols refer to simulation data. The same sequence of c/> 
values is shown in both panels. 



observed between the two studied isochores (j) = 0.288 
and = 0.314 . While in the (j) = 0.288 case P(T) de- 
creases continuously on cooling, in the (j) — 0.314 case 
the low T behavior of P is reverses and P approaches 
a positive finite value on cooling. This different low-T 
trends indicated that for < 0.3, on cooling the network 
becomes stretched (negative pressures), in the attempt 
to preserve the connected bonded state. This implies 
that at low T, there is a driving force for phase sepa- 
rating into a fully connected unstressed network and a 
gas phase. This also suggests that the spinodal curve 
ends at T = around (p = 0.3. At (/> « 0.3, the pack- 
ing fraction is optimal for the formation of an unstressed 
fully connected network at low T. The bond formation 
on cooling does not require any stretching and it reverses 
the T-dependence of P. (iii) Between 0.3 ^ 4> ^ 0.38 a 
minimum of P appears. The existence of a minimum in 
P{T) along isochores evidences the presence of density 
anomalies (i.e. expansion on cooling along isobars) since 
points in which {dP/dT)v = 0, by a Maxwell relation, 
coincide with points in which a = {dV/dT)p = 0, i.e. 
with points in which density anomalies are present. The 
simplicity of the model allows us to access the different 
contributions to P and investigate the origin of the in- 
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FIG. 6: Components of the pressure a.t (p — 0.314. The total 
P is decomposed in ideal gas, HS and bonding components. 
Note the isochoric minimum in P around T — 0.105, a signa- 
ture of an isobaric density maximum. 

crease of P on cooling. In the PMW, apart from the 
trivial kinetic component contribution, the only positive 
component to P arises from the HS interaction. Inter- 
estingly enough, the HS component increases on cooling. 
Such an increase in the HS repulsion, indirectly induced 
by the formation of the bonding pattern, in the range 
0.30 ^ (j) ^ 0.36 appears to be able to compensate the 
decrease in the bonding component of P. 

To confirm the presence of density anomalies it is 
instructive to look at the V dependence of P along 
isotherms, shown in Fig. [7| Again, the simulation data 
are consistent with the Wertheim theory predictions only 
at large T and indeed it was already noted that no density 
anomalies are found within the theoryP|. The simulation 
data also show a clear crossing of the isotherms around 
a volume per particle v = lA and 1.7, corresponding to 
(j> = 0.314 and 4> = 0.38. Again crossing is indicative of 
the presence of density anomalies. The increase of P on 
cooling, between (/> — 0.314 and (j) — 0.38 suggest also a 
possible emergence of a second Van der Waal-type loop 
(in addition to the gas-liquid one) for T lower than the 
one we are currently able to equilibrate. The possibil- 
ity of a second critical point between two liquid phases 
of different densities has been discussed at length in the 
past [33, following the discovery of it|^ in one of the 
first models for water [39j. 



C. g{r) 

The PMW radial distribution functions for T > 0.15, 
have been reported previously ^Ij. Here we focus on the 
interesting structural changes observed during the devel- 
opment of the bond network in goo and gu-LP, a T- 
region which was not possible to access in the previous 
simulations. The 1700 provides information on the cen- 
ter to center particle correlation while gn-Lpif) contains 




b 1 I 1 I 1 I 1 I 1 d 

1 1,5 2 2,5 3 3,5 




FIG. 7: Isotherms of P according to the Wertheim theory 
(top) and as calculated from the simulation data (bottom) 
as a function of the volume per particle v = . Symbols 
refer to simulation data. Note the crossing of the different 
isotherms ai v — 1.4 and 1.7. 



information on the bonding and on the attractive com- 
ponent of the pressure. 

Fig. IHl shows goo{i') at three different packing frac- 
tions. In the interval 1 < r < 1.1 the function is highly 
peaked, a consequence of the distance imposed by bond- 
ing. Outside the bonding distance (r > 1.1), goo{r) 
shows significant oscillations only at low T. A peak, be- 
side the bonding one, is observed at r « 1.7 correspond- 
ing to the characteristic distance between two particles 
bounded to the same central particle in a tetrahedral 
geometry. The absence of the information about the ge- 
ometry of the bonding sites in the theory of Wertheim 
is responsible for the absence of the peak at 1.7tT and 
the breakdown of the predictive ability of the Wertheim 
theory as soon as a particle is engaged in more than two 
bonds. A few observations are in order when comparing 
the (j) dependence of gooif)- At low 0, the tetrahedral 
peak at r « 1.7 is the only peak in goo{T)- When (p ap- 
proaches the optimal network density a clear tetrahedral 
pattern develops and goo{f = 1-7) becomes larger than 
two. The tetrahedral peaks at « 1.7 is followed by oscil- 
lations extending up to 4(7. At even larger 0, there is still 
a residual signature of tetrahedral bonding at 1.7(7, but 
the depletion region for r > 1.1(7 is not developed any 
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FIG. 8: Particle-particle radial distribution function at 
0.105, 6 = 0.288, S = 0.380. 



longer, signaling a competition between the HS packing 
(which favor a peaks at positions multiple of cr) and the 
local low density required by bonding. 

Fig. H compares, at = 0.314, the 00, HH and H- 
LP radial distribution functions in linear scale. In all 
three functions, the progressive structuring induced by 
the bonding is clearly evident. Even gHnir) shows very 
clear signs of spatial correlations, which are induced by 
the tetrahedral geometry of the bonding and by the ge- 
ometry by which the bonding between H and LP prop- 
agates. Indeed, in the PMW model the interaction be- 
tween different H sites is zero. 
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FIG. 9: Radial distribution functions for OO, HH and H — 
LP pairs at the optimal network density — 0.314. Insets in 
goo and Qh-lp provide enlargements of the contact region. 
On cooling a significant structure appears, associated to the 
intense bonding. 



D. S{q) 

The structure factor of the system, defined in term of 
the particle's center coordinates ri as, 



1 



N 



i=l 



(10) 
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provides information on the wave vector dependence of 
the density fluctuations. In isotropic systems, S{q) is 
function of the modulus q. The behavior of S{q) at small 
q provides indication on the phase behavior, since an in- 
crease of S{q) at small q indicates the development of 
inhomogeneities with length-scale comparable to the sys- 
tem size studied. As an indicator of the location of the 
phase boundaries (of the liquid-gas spinodal line) , we es- 
timate the locus of points in (T, 0) where S{q) for the par- 
ticles centers becomes larger than 5 at small q. This locus 
is reported in Fig. 0] For cj) > 0.28 S{q) does not show 
any sign of growth at small q in the region of T where 
equilibration is feasible, being characterized by values of 
S{q) at small q of the order of 0.1. This confirms that, at 
this packing fraction, there is no driving force for phasi 
separation, since the average density has reached a valui 
such that the formation of a fully connected network o 
bonds does not require a local increase of the packin; 
fraction. It is also important to stress that dX (j) — 0.288 
at the lowest studied T, the average number of bond pe: 
particle is 3.8, and hence the system is rather close to iti 
ground state and no more significant structural change; 
are expected on further cooling. 

Fig. Uni shows S{q) at </) = 0.105, <p = 0.288 anc 
(f) — 0.385. The = 0.105 case has been chosen to shov 
the significant increase in S{q) associated to the approacl 
of the spinodal curve. The case = 0.288 shows botl 
the absence of a small g-vector divergence and the clea: 
development of the typical g-pattern of tetrahedral net 
works. On cooling the peak at qa — 27r, characteristic o 
excluded volume interactions splits in two parts. A pre 
peak around gcr « 5 and an intense peak around gcr w 8 
The case = 0.385 confirms that the packing fraction i: 
now so high that a full tetrahedral network cannot de- 
velop and the splitting of the main peak in two distinct 
components is very weak and visible only at the slowest 
investigated T. 



E. Percolation 

The PMW, as all other models based on HS and SW in- 
teractions, is particularly suited for calculation of bond 
properties, since a bond between particle i and j can 
be unambiguously defined when the pair interaction en- 
ergy between i and j is —uq. In the case of continuous 
potentials such a clear cut bond definition is not possi- 
ble and several alternative propositions have been put 
forward liol, IS']. We focus here on the connectivity prop- 
erties of the equilibrium configurations. We use standard 
algorithms to partition particles into clusters. Configura- 
tions are considered percolating when, accounting for pe- 
riodic boundary conditions, an infinite cluster is present. 
More explicitly, to test for percolation, the simulation 
box is duplicated in all directions and the ability of the 
largest cluster to span the replicated system is controlled. 
If the cluster in the simulation box does not connect with 
its copy in the duplicated system then the configuration 
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FIG. 10: Particle-particle structure factor &t (j) = 0.105, cj> = 
0.288, (j) = 0.3858. Note that &t 4> ^ 0.105, an intense signal 
develops at small q, related to the approach to the spinodal 
instability. Small q intensity is completely missing at the 
higher (j) shown. 



is assumed to be non-percolating. The boundary between 
a percolating and a non-percolating state point has been 
defined by the probability of observing infinite clusters 
in 50% of the configurations. The resulting percolation 
line is reported in Fig. ^ State points on the right side 
of the line are characterized by the presence of an infinite 



9 



cluster. Still, at this level of definition, percolation is a 
geometric measure and it does not provide any informa- 
tion on the lifetime of the percolating cluster. 

The percolation line, like in simple SW potentials, 
crosses the spinodal curve very close to the critical point. 
Differently from the SW, the percolation locus does not 
extend to infinite T, since at high T, even at large 0, the 
reduce particle surface available for bonding prevents the 
possibility of forming a spanning network with a random 
distribution of particles orientations. Along the perco- 
lation line, about 1.5 bonds per particle are observed, 
with a small trend towards an increase of this number on 
decreasing (p. In terms of bond probability pf,, this cor- 
respond to Pi, « 0.375, not too different from the bond 
percolation value of the diamond lattice, known to be 
0.388|42j. 



IV. DYNAMICS 

Thermodynamic and static properties of the PMW 
presented in the previous section clarify the location of 
the regions in which the bond network forms, the region 
where the liquid-gas phase separation takes place and the 
region at high </> where packing phenomena start to be 
dominant. In the following we present a study of the 
diffusion properties of the model in the phase diagram, 
with the aim of locating the limit of stability of the liquid 
state imposed by kinetic (as opposed to thermodynamic) 
constraints. 



A. MSD 

We focus on the mean square displacement < r^{t) > 
of the particle centers, as a function of T and (p, calcu- 
lated from the newtonian dynamic trajectories. Fig. ITTI 
shows < r^it) > for a few selected isochores. For short 
time < r^{t) >=< > , where < v"^ >= 3/2kBT is 
the thermal velocity. At high T, the short-time ballistic 
behavior crosses over to a diffusion process (< t) 
directly. At low T, the ballistic short-time and the diffu- 
sive long time laws are separated by an intermediate time 
window in which < r^(i) > is approximatively constant, 
an indication of particle caging. 

Several features of < r^{t) > are worth pointing: (i) 
For <j) < 0.209, the spinodals are encountered on coohng 
before the caging process is visible. The phase separa- 
tion process sets in well before particles start to feel the 
caging process, (ii) The static percolation curve reported 
in Fig. ^ has no effect on dynamics. There is no dynamic 
arrest at the static percolation transition, (iii) For 
such that a well developed tetrahedral network can form, 
it is possible to cool the system down to temperatures at 
which, on the scale of simulation, arrest is observed, in 
the absence of any phase separation. < r^{t) > develops 
a clear intermediate region where only the dynamic inside 
the cage is left. At this p, the caging is not associated to 



T = 0.12 
A^T = 0.13 
<-<iT = 0.14 

T = 0.17 
-1—1-1 = 0.20 




10"^ 10-' 10° 10' 10^ 10^ 10"* 10^ 
t 

10^1 





10' 




10" 










V 


10-' 



G-G T=0.095 

[3-e:t=o.io 

T=0.11 
A-^T=0.115 

T=0.12 

T=0.13 
>-t>T=0.135 
H — hT=0.14 

T=0.15 
*-*T=0.17 

T=0.20 



TTi — I I 1 1 iiiii — I I 1 1 mil — I I 1 1 Mill — jm 1 1 III 




10- 



10-^^ 

10-2 10' 10" 10' 10^ 10^ W 10^ 
t 

FIG. 11: Mean square displacement for different Ts, at three 
different (j> values. Top: (p = 0.131, Center: <f) — 0.314, Bot- 
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excluded volume interactions, but to the formation of en- 
ergetic bonds 4§\ . (iv) The plateau value in < (t) > is 
a measure of the localization length induced by the cage. 
To visualize the dependence of the localization length, 
we show in Fig. ^] < {t) > for three different state 
points (0-r) with the same long time diffusivity. The 
cage length is always significantly larger than the typical 
HS value (< r^(i) >~ 0.01) and grows on decreasing p. 
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FIG. 12: Mean square displacement along a constant D path. 
Note the (j) dependence of the plateau at intermediate times, 
which provides an estimate of the caging length. 



B. Diffusion Coefficient 

The long time limit of < r^(t) > is, by definition, 6Dt, 
where D is the diffusion coefficient. The 4> and T depen- 
dence of D is shown in Fig.^| We show log(I?) both vs. 
T and vs 1/T. Again, a few considerations are in order: 
(i) The range of D data covers about five orders of magni- 
tude. The data for (j> < 0.24 are limited in T by the phase 
separation process, while the data for cj) > 0.26 are lim- 
ited by computational resources, since equilibration can 
not be reached within several months of calculations, (ii) 
Data for cf) > 0.26 crosses around T « 0.105, suggesting 
a non monotonic behavior of the (j) dependence of the dy- 
namics, (ii) The early decay oiD with T can be described 
with a power-law \T — TmctP- Power law fits, limited 
to the region of T between T = 0.11 and T = 0.15, cover 
the first two-three orders of magnitude in D, in agree- 
ment with previous studies of more detailed models for 
water |43l l44l l45l | and with the previously proposed MCT 
interpretation of theml4^ |4^ |47|, li^- (iii) A cross-over 
to an Arrhenius activated dynamics is observed at low 
T. Activated processes become dominant in controlling 
the slowing down of the dynamics. The activation en- 
ergy is « 4 close to the optimal network (j), suggesting 
that at low T diffusion requires breaking of four bonds. 
The cross-over from an apparent power-law dependence 
to an Arrhenius dependence has also been observed in 
simulations of other network forming liquids, including 
silica^^ IH^I and more recently water [5l|. The low T 
Arrhenius dependence also suggests that in the region 
where bonding is responsible for caging the vanishing D 
locus coincides with the T = line. 

Particularly interesting is the behavior of D{(p) along 
isotherms. An almost linear dependence at small (f> (up 
to (j) = 0.235) is followed by a non monotonic behav- 
ior. Below T — 0.11, a diffusion anomaly is observed 
in the T and 4> region where the tetrahedral network 
develops. Around (j) = 0.34 an isothermal compression 




FIG. 13: Temperature dependence of the diffusion coefficient 
along isochores. The dashed line is an Arrhenius dependence 
with activation energy equal to 4. 



of the system generates a speed up of the dynamics. 
Above (j> w 0.35, D starts to decrease again on increas- 
ing packing. Diffusivity anomalies of the type observed 
in the PMW are found in several tetrahedral network 
forming liquids, including water [s^. The explanation for 
this counterintuitive (jj dependence of the dynamics is to 
be found in the geometric constraints requested by the 
tetrahedral bonding requiring an open local structure. 
Increasing (j) destroys the local bonding order with a re- 
sulting speed up of the dynamics. 



C. Isodiffusivity (and arrest) lines 

A global view of the dynamics in the {T~(j)) plane is of- 
fered by the isochronic lines, i.e. the locus of state points 
with the same characteristic time|53|. In the present 
case we focus on the isodiffusivity lines. The shape of 
the isodiffusivity lines, extrapolated to I? ^ provides 
a useful indication of the shape of the glass transition 
line 13 EH 113 • Fig. [T31 shows the isodiffusivity lines for 
several different values of D, separated each other by one 
order of magnitude. The slowest isodiffusivity lines are 
only weakly T dependent at low 4>. For small values of 
D, iso-diffusivity lines start from the right side of the 
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FIG. 14: Diffusion coefficient along isotherms. Note the non- 
monotonic behavior which develops for T < 0.11. 



spinodal, confirming that slow dynamics is only possible 
for states with <j) > (j>c- At large the isodiffusivity lines 
bend and become parallel to the T axis, signaling the 
cross-over to the hard-sphere case. Extrapolating to zero 
the T (or cj)) -dependence of D it is possible to provide es- 
timates of the dynamic arrest line. In the present model, 
the low T-dependence of D along isochores is well mod- 
eled by the Arrhcnius law and hence technically arrest is 
expected at T = 0. The shape of the iso-diffusivity lines 
suggests that the vertical repulsive glass line (controlled 
by excluded volume effects) starting at high T from the 
HS glass packing fraction meets at a well defined the 
T = bond glass line. 

The shape of the PMW isodiffusivity lines is very sim- 
ilar to the short-range square well case, for which a flat 
T-independent "attractive" glass line crosses (discontin- 
uously) into a perpendicular (j) independent "repulsive" 
glass line|55l f57l|. Differently from the SW case, in the 
PMW the equivalent of the attractive glass line extends 
to much smaller <j) values, since the reduced valency has 
effectively reduced the space in which phase separation 
is observed^. It is also worth pointing that the shape of 
the isodiffusivity lines at low cf) is similar to the shape 
of the percolation line. As in all previously studied 
models 7, 55], crossing the percolation line does not co- 
incide with dynamics arrest, since the bond lifetime is 
sufficiently short that each particle is able to break and 
reform its bonds. 



D. D vs. E - Egs 

At the optimal network density, the low T behavior 
of both D and E — Egg (which, as discussed above, is 
also a measure of the number of broken bonds) is Arrhe- 
nius. This suggests to look more carefully in the relation 
between the activation energy of the two processes. One 
possibility is offered by a parametric plot of D vs E — Egs 
in log-log scale, so that the slope of the straight line pro- 



FIG. 15: Isodiffusivity lines in the (T — (p) plane. An excur- 
sion of five orders of magnitude in D values is explored. All 
lines start from the spinodal and end at infinite T at the cor- 
responding HS location. At small D, lines cannot be contin- 
ued above (p = 0.5 since there the HS interaction is dominant 
and the system crystallizes. Extrapolating along isochores the 
observed Arrhenius functional form suggest an ideal D = 
arrest line at T = 0. 



vides the ratio of the two activation energies. Such a plot 
is shown in Fig. ^| We find the remarkable results that 
close to the optimal network (j>, the slope of the curve has 
exponent four, i.e. D ^ (1 — pt)^, where pi, is the prob- 
ability that one of the four bonds is formed (and hence 
l—pb is the probability that one of the four possible bonds 
is broken), suggesting that the elementary diffusive pro- 
cess requires the breaking of four bonds. A functional 
law for diffusion in a tetrahedral model of this type was 
proposed by Teixera(58j to interpret the T dependence 
of D in water in the context of the percolation model 
developed in Ref. [s^- A similar dependence has been 
recently reported for a model of gel-forming four-armed 
DNA dendrimers[60|. 




FIG. 16: Diffusion Coefficient vs E—Egs for different (j> values. 
The dashed line is a power-law with exponent four. 
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E. D - MD vs. MC 

All dynamic data presented above refer to event-driven 
newtonian dynamics. Indeed, Monte Carlo simulations 
intrinsically miss dynamic informations being based, in 
their simpler formulations, on random displacements of 
the individual particles. Still, if the random displacement 
in the trial move is small as compared to the particle 
size the sequence of MC steps can be considered a pos- 
sible trajectory in configuration space. When this is the 
case, the number of MC-steps (each step being defined 
as an attempted move per each particle) plays the role 
of time in the evolution of the configurations in config- 
uration space. In the absence of interactions, a particle 
evolved according to the MC scheme diffuses with a bare 
diffusion coefficient D^jfj fixed by the variance S^^j of 
the chosen random displacement along each direction (in 
our calculations we have used an uniform distribution of 
displacements with a variance of i5^(^ = (0.1)^/12, cor- 
responding to -D°/c — 3(5|/(7/6 in units of u^/MC-step). 
If needed, -D^v/c provides a mean to associate a physi- 
cal time to the MC-step. At low T, when slow dynamic 
processes set in (favored by bonding or by packing), it 
is expected that the microscopic dynamics becomes irrel- 
evant (except for a trivial scaling of time). The escape 
from the cage created by neighboring particles is indeed 
a much rare event as compared to the rattling of the par- 
ticles in the cage. Under these conditions, the slow dy- 
namic processes become independent on the microscopic 
dynamics, and hence Newtonian, Brownian and MC show 
the same trends. FigEI shows that this is the case for 
three (f> values. In all cases, at low T, the T dependence 
of D^^'-^ and D is identical. Moreover, the scaling fac- 
tor between MC and MD dynamics is independent of 
4>, suggesting that at low T, with the chosen units, the 
relation D^'^ / Dl^^j = holds. From comparing MC 
and MD data we find that the proportionality constant 
^ « 10 and shows no state-point dependence. To con- 
firm that caging is fundamental to observe independence 
of the slow-dynamics from the microscopic one, we look 
at the shape of < r^(t) > (Fig. [TT|l . finding that at the 
T at which MC and MD dynamics start to coincide a 
significant caging is present. 

Since the microscopic time of the MC dynamics is not 
affected by temperature (being always fixed by the vari- 
ance of the random displacements) it is interesting to 
consider the relation between D and E — Egs also for 
D^'-' , shown in Fig. ^| at the optimal network density 
(j) = (f) = 0.314. Again, the slope of the curve has ex- 
ponent four, but compared to the MD case, the region 
of validity of the power-law covers the entire range of T 
studied, from very high T (where the number of bonds is 
negligible) down to the lowest equilibrated temperature, 
covering more than 4 order of magnitude. The valid- 
ity of the relation D ~ (1 — p^Y extends up to high T, 
when the system is well above percolation and there is 
no evidence of a tetrahedral network (as shown in the 
structural data reported in Fig. 1101 andl5|). The extended 
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FIG. 17: Comparison between the MD and MC diffusion 
coefficient at three diflferent 4> values. The MC data are also 
shown multiplied (by a common factor 0.1) to better visualize 
the low T overlap. 



validity of the power-law, with an exponent exactly equal 
to the valence of the model is highly suggestive and, in 
principle, very important for theoretical considerations, 
since it appears to cover either the region of tempera- 
ture where liquid dynamics is observed, either the low T 
states where signatures of slow dynamics (see Fig lll|l are 
very well developed. The limit of validity of this finding 
needs to be carefully checked in other primitive models 
with different valence and with more realistic models of 
network forming liquids. 
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FIG. 18: Relation between D , normalized by the bare 
MD diffusion constant -Djv/c ''■nd E — Egs for MC dynamics. 
Note that the MC data follow over more than five orders of 
magnitude a simple fourth-power law (full red line). 



V. CONCLUSIONS 

Results presented in this manuscript cover several ap- 
parently distinct fields. 

To start with, results presented here can be discussed 
in relation to the dynamic and thermodynamic proper- 
ties of water. We have shown that the thermodynamic of 
the PMW includes, beside the compressibility anomalies 
reported beforeTl, also density anomalies (at much lower 
T). The source of the density anomalies is shown to be 
associated to the establishment of the bond network in 
the tetrahedral geometry. On cooling (along isochores) 
the energetic driving force which favors the formation 
of the bond, due to geometric constraints associated to 
the formation of the open tetrahedral structure, forces 
the pressure to increase and hence generating a density 
maximum state point. The simplicity of the PMW al- 
lows us also to clearly detect an optimal network density, 
at which the ground state of the system (i.e. the state 
in which each particle is involved in four bonds) can be 
closely approached. At this (/> the T-dependence of the 
potential energy is the most pronounced, generating a 
minimum in the isothermal (j) dependence. The pres- 
ence of a minimum in E[(j))\T is highly suggestive since it 
indicates the possibility of a liquid- liquid phase sep- 
aration at T lower than the one we have been able to 
equilibrate. We have also shown that at this optimal 0, 
low T dynamics slows down with the fourth power of the 
probability of broken bonds, i.e. the dominant compo- 
nent to dynamics arises from single particle motions, and 
specifically of the particles which happen to have all four 
bonds broken at the same time. We have also shown that, 
like in real water, diffusion anomalies are observed. At 
low T, the decrease of the diffusivity on increasing (j) is re- 
versed once the optimal network density is reached. For 
higher cf), the progressive destruction of the bond network 
due to the increased packing fastens the dynamics. For 



even higher (j), D((j)) resumes its standard decreasing be- 
havior associated to the approach of the excluded volume 
glass transition. Diffusion and density anomalies in the 
PMW models are thus strongly related, similarly to what 
has been observed in more realistic models for water [6^. 
The simplicity of the model is crucial in clarifying these 
aspects since the hard-core and square well interactions 
guarantee the absence of volumetric effects related to the 
T-dependence of the vibrational amplitudes. 

A second interesting aspect of the presented results 
concerns the dynamics in network forming systems. The 
present study provides a complete characterization of the 
dynamics in the entire {(p — T) plane, from the smallest 
possible (j) liquid state points up to the close packed state. 
From the reported data, the relative role of the energy 
and of the packing in controlling the dynamics stands up 
clearly. The isodiffusivity lines are essentially parallel to 
the (/)-axis (i.e. T controlled) in the network low cf) region 
and are essentially parallel to the T-axis (i.e. 4> con- 
trolled) at larger (f>. Interesting enough, along isochores, 
low T dynamics follows and Arrhcnius law, the landmark 
of strong-glass forming behavior 63, 64J. The Arrhenius 
law is foreseen by a T region where dynamics has a strong 
T dependence, compatible with a power-law dependence. 
In this power-law region the first signatures of caging 
in the mean square displacement are observed. Similar 
changes in the dynamics have been observed in previ- 
ous studies of siHca^^ [s^, 113 , water [sij and sihcon[6^. 
In particular, for the case of silica and water, it has been 
suggested that the region where dynamics start to feel the 
presence of energetic cages can be interpreted in terms of 
mode coupling theoryjH III lei IgI IgI IgI Fol l7l| . 

Dynamics at the optimal network (j) is particularly sug- 
gestive. Although in the present model, slowing down of 
the dynamics prevents equilibration of the supercooled 
liquid to very low T, at the lowest T simulations the av- 
erage number of bonds has gone up to 3.8 per particle. 
In this respect, further structural and dynamic changes 
are hard to foresee. This suggests that the Arrhenius 
behavior is retained down to T = 0. Such speculation 
is reinforced by the numerical values of the activation 
energy of D which is found to be « 4uo, i.e. correspond- 
ing to the breaking of four bonds. This suggests that 
in network liquids, the limited valency imposed by the 
directional forces fixes a well defined energy of the local 
configuration and a discrete change of it which is reflected 
in the Arrhenius behavior. The presence of a limited va- 
lency and a well defined bond energy scale appears to 
be the key ingredient of the strong liquids behavior*^. 
It is also worth to explore in future works the possibil- 
ity that the optimal network density plays, in studies of 
one component systems, the same role as the reversibility 
window 72] in bulk alloy glasses. Connections with the 
concept of self-organization in network glasses [t^ should 
also be pursued. 

A further aspect of this work concerns the relative lo- 
cation between the liquid-gas spinodal and the kinetic ar- 
rest lines, whose shape is inferred by the study of the isod- 
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ifFusivity lines. As in the short range SW modellzllzi, 

the kinetic arrest lines ends in the right side of the spin- 
odal, i.e. in the liquid phase. But differently from the 
SW case, the limited valency has shifted the right side of 
the spinodal to very small (/) values, <j> < 0.25. Indeed, the 
limited valency effectively disfavors condensation of the 
liquid phase reducing the driving force for phase separa- 
tion and making it possible to generate low packing frac- 
tion arrested states in the absence of phase separation, 
i.e homogeneous single phase stable in equilibrium, at 
low r[7^. The possibility to access low T homogeneous 
supercooled states for (j> > 0.25 characterized by a glassy 
dynamics, driven by the bonding energy as opposed to 
packing, confirms the findings of the zero-th order model 
with limited valency reported in Ref. ^ . The absence of 
geometric correlation between the bonding sites, the key 
ingredient of the maximum valency model 7] is thus not 
crucial for the stabilization of the network. The role of 
the geometric constraint appears to be the reduction in 
the range of (j) values where the fully bonded disordered 
state can be reached. Two different arrest mechanisms 
characterize the dynamics of network systems. Arrest 
due to the formation of energetic cages, with an arrest 
line which runs almost parallel to the (p axis, and arrest 
due to excluded volume effects, with an arrest line par- 
allel to the T axis. These two lines are reminiscent of 
the attractive and repulsive glas s lines observed in short- 
range attractive colloids Issl IstL l6^ l77| . Connecting 
the results presented in this article with previous studies 
of network forming liquids^JS, ,49J, it is tempting to spec- 
ulate that mode-coupling theory predicts satisfactory the 
shape in the {(j) ~ T) plane of the dynamic arrest lines. 
Still, while in the region where excluded volume controls 
caging the relative error in the location of the glass line 
is limited , in the case in which bonding mechanism is 
dominant in generating arrest, the location of the MCT 
line can be significantly distant from the actual dynamic 
arrest line (technically located at T = 0, being dynamics 
Arrhenius), due to the role of activated bond-breaking 
processes which offer a faster channel for the decay of 
the correlations. The evaluation of the MCT lines for the 
present model, in principle feasible within the site-site ap- 
proach developed by Chong and Goetze^,^7^ or within 
the molecular approach developed by Schilling | 46ll80| can 
help clarifying this issue. 

The possibility of an intersection between the excluded 
volume arrest-line (starting at high T from the HS glass 
packing fraction ) and the bond-controlled T = ar- 
rested line is particularly suggestive. The shape of the 
iso-diffusivity lines supports the possibility that the ver- 
tical repulsive glass line meets at a well defined (j) the 
T = bond-controlled glass line. If this scenario is cor- 
rect and general, one would conclude that the fragile and 
strong kinetic behavior is intimately connected to the 
dominant mechanism of arrest (fragile for excluded vol- 
ume and strong for bonding) and, more interestingly, that 
strong behavior can only be observed when the interac- 
tion potential is such that less than six neighbors are 



present (i.e. in network forming systems). Indeed, only 
under these circumstances the suppression of the liquid- 
gas phase separation makes it possible to approach the 
T ~ bond-controlled glass line. 

An additional comment concerns the relation between 
gel and glass arrest states. Results reported in this ar- 
ticle confirm, one more, that in this class of models 
the geometric percolation line does not have any im- 
pact on the dynamic arrest, since at percolation the life- 
time of the bond is still rather small. Only when the 
system is well inside the percolation region, the bond 
lifetime has slowed down significantly to affect all mea- 
surements of global connectivity with an intrinsic time 
scale shorter than the bond lifetime (as for example fi- 
nite frequency shear viscosity). Indeed, already long time 
ago it was noted for the case of water's^ that bond 
percolation is irrelevant to any thermodynamic or dy- 
namic anomaly. More sophisticated models, incorporat- 
ing bond-cooperativity or significant entropy contribu- 
tions to bonding (as the case of polymeric gels) may re- 
duce the differences between dynamic arrest states and 
percolation [60i|. 

Despite the difference between percolation and arrest 
lines, if one consider the present model as a system of 
colloidal particles with sticky interactions, one would be 
led to call the arrested state at 0.3 ^ (f> ^ 0.5 a gel, led 
by the fact that the arrested state has a low (j) open con- 
nected structure. Similarly, if one consider the PMW as a 
model for a network liquid, one would be led to name the 
same arrested state a network glass. While we cannot of- 
fer any resolution to this paradox with the present set of 
data, future work focusing on the shape of the wavevec- 
tor dependence correlation functions and the resulting 
non ergodicity parameters can help clarifying this issue 
and confirm/dispute the hypothesis on the differences 
between gels and glasses recently proposed^ [t^, Isij . 
At the present time, we can only call attention on the 
fact that a continuous change from energetic cages to 
excluded volume cages takes place on increasing (f). 

A final comment refers to the propensity of the sys- 
tem to form disordered arrested states. Despite the rel- 
evant amount of supercooling^^, in all studied state 
points where a network structure is present, we have 
not observed any sign of crystallization. The kinetic 
suppression of the crystallization phenomenon can be 
traced to the similar energy characterizing the crystal 
state and the fully bonded disordered state, vanishing 
the energetic driving force toward crystallization. The 
observed propensity to form gel states as opposed to crys- 
talline states manifested by the studied model (which 
can be seen also as a model for short-range sticky col- 
loid particles as well as globular proteins with aeolotopic 
interactions 0) may well explain the difhculty of crystal- 
lizing some class of proteins. It also warn us about the 
relevance of the dynamic arrest phenomenon in the future 
attempts to build a colloidal diamond photonic crystal, 
made of particles with short-ranged patchy interactions. 
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VII. APPENDIX: AN EVENT-DRIVEN 
ALGORITHM FOR HARD SPHERES WITH 
PATCHES. 

In an event driven (ED) algorithm, events such as 
times of collisions between particles and cell crossing have 
to be taken into account. All these events have to be or- 
dered. Code must be written in such a way that locating 
the next event and insertion/deletion of new events have 
to be performed efficiently. In literature, several ED al- 
gorithms for simulating hard-sphere systems exist and 
several propositions on how to handle such events effi- 
ciently have been reported. One elegant approach, pro- 
posed twenty years ago by Rapaport :82| , arranges events 
into an ordered binary tree (calendar of events) so that 
insertion, deletion and retrieving of events can be done 
with an efficiency O(logiV), 0(1) and O(logA^) respec- 
tively, where N is the number of events in the calendar. 
We adopted this solution to handle the events calendar in 
our simulation, adding only a redefinition of event time in 
order to avoid round-off problems which are found when 
extremely long simulation runs are performed. 



A. Motion of rigid bodies 

The orientation of a rigid body can be conveniently 
represented by the 3 column eigenvectors (with i = 
1,2,3) of the inertia tensor expressed in the laboratory 
reference system. These vectors form an orthogonal set 
and can be arranged in a matrix R, i.e. 



R = *(uo ui U2) 



(11) 



where *A indicates the transpose of the matrix A. This 
matrix is such that if x are the coordinates of the labo- 
ratory reference system and x' are the coordinates of the 
rigid body reference system, it turns out that: 



Rx 



(12) 



In what follows, we assume that the three eigenvalues 
of the inertia tensor are all equal to /. Naming w — 
{wx, Wy,Wz) the angular velocity of a free rigid body, the 
matrix 17 is defined as 



(13) 





( 




Wy 


n = j 







—w._ 












Knowing the orientation at time t 
R(<) at time t is: t83i>.84|: 

Kit) =R(0)(I + M) 



0, the orientation 



where M is the following matrix: 
sin{wt) _ 1 



M: 



-n 



cosjwt) ^2 



w 



(15) 



and w = ||w||. Note that if w = then R(i) = R(0). To 
derive Eq. ^1 consider that: 

'Rit) = (Ui(t) U2(0 U3(t)) (16) 

= (*(I + M)ui *(I + M)u2 *(I + M)U3 ) 

where we remember that are column vectors. Hence, 
if w = wh, we have after some algebra: 

Ui{t) = U; • n n -f cos(wi)(ui — n • n) + sm{wt) n x 

(17) 

that is the so-called Rodriguez's formula or Rotation for- 
mula, i.e. a rotation of an angle wt around the axis n. 
To conclude if one has to update position and orienta- 
tion of a rigid body, that is freely moving, this can be 
accomplished doing: 



x(i) = x(0) -t- vt 



nit) = R(0)(I + M) 



(18a) 



(18b) 



where x(<) is the position of the center of mass of the 
rigid body at time t and v is its velocity. 



B. Hard-Sphere with interacting patches 

In the present model, each particle is modeled as an 
hard sphere with n spherical patches arranged in fixed 
site locations. In the present case, the site-site interaction 
is a SW potential. 



usw 



-Uq 







ii r < S 
otherwise 



(19) 



where 6 and uq are the width and the depth of the SW 
. For the following discussion, the SW interaction can 
be visualized as a sphere of diameter S centered on the 
site location. Similarly, one can visualize the particle as 
a rigid body composed by the hard-sphere joined to the 
spheres located on the sites. In what follows, we identify 
a particle with the resulting surface. Defining distance 
dAB between two particles A and B as the shortest line 
connecting two points on distinct particles, i.e. 



dA 



min di 



iA,iB 



(20) 



where iat^b € {0, ...n} and labels the hard sphere, 
1 . . . n label the n spherical patches and dij^i^ is the dis- 



(14) tance between the two spherical patches iA and iB- 
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C. Prediction of time-of-collision 

1. Finding the contact time 

We separate the collisions between two particles in the 
hard-sphere part of the potential and the site-site inter- 
action part. The time of collision t^s between the hard 
sphere cores can be evaluated as usual [s^l ■ The smallest 
time of collision among all spherical patch pairs is tgt ■ 
Time of collision of the two particles is 



tc = inm{ths,tst} 



(21) 



To find the time-of-collision of two interacting patches we 
assume that it is possible to bracket it. I.e., we assume 
(see further subsections) that the time of collision tgt is 
such that ti < tst < t2 where the product d{ti)d{t2) < 0. 
Thus, the "exact " time of collision is provided by the root 
of the following equation: 



\hAt)-r^At)\\^s 



(22) 



where r^^ and r.^^ are the two site locations. 



2. Linked lists and Centroids 

As described in [131 to speed up a ED molecular dy- 
namics of hard spheres one can use linked lists. For 
a system of N identical particles inside a cubic box of 
edge L, we define the "centroid" [s^ [s^ as the smallest 
sphere that contains the particle (the HS and the spheri- 
cal patches). Linked lists of centroids may be quite useful 
to reduce the number of objects to check for possible col- 
lisions and in addition they can be used to restrict the 
time interval within which searching for the collision. We 
divide the box into cells so that into each cell con- 
tains at most one centroid. After that we build the linked 
lists of these centroids and handle these lists as done in a 
usual ED molecular dynamics of hard spheres [s^ • This 
means that whenever an object cross a cell boundary one 
has to remove such an object from the cell the particle 
is coming from and add this object to the cell where it's 
going to. 

Now consider that one has to predict all the possible 
collisions of a given particle, which is inside a certain cell 
m. As for the hard spheres case we take into account only 
the particles inside the adjacent cells (see for more 
details) and we predict the times of collisions with these 
objects. Consider now two particles A and B at time 
t = and their centroids Ca and Cb- Three possible 
cases arise: 

1. Ca and Cb do not overlap and, evaluating their 
trajectory no collision between the two centroids is 
predicted. In this case A and B won't collide as 
well. 

2. Ca and Cs do not overlap but they will collide: in 
this case, we calculate two times ti and t2, brack- 
eting the possible collision between A and B: ti is 



defined as the time when the two centroids collide 
and start overlapping and t2 is the time when the 
two spheres have completely crossed each other and 
do not overlap any longer. 

3. Ca and Cb overlap: in this case ti = and t2 is 
defined as the time at which the two centroids stop 
overlapping. 



3. Fine temporal bracketing of the contact time 

Here we show how a refined bracketing of solution of 
Eq. |211can be accomplished. First of all we give an over- 
estimate of the rate of variation of the distance between 
two patches za and Ib, that is: 



< 
< 



d 

dt 



r 



< 



W-tAiB I 

\Yab+(^a X (Vi^ 



Ra) - t^s X (r; 

\\u;B\\LB = dZZ 



^b)\\ 
(23) 



where the dot indicates the derivation with respect to 
time, , Tig are the positions of the two sites with re- 
spect to laboratory reference system, v^^^^, is the rela- 
tive velocity of the two sites, V^b is the relative velocity 
between the centers of mass of the two particles, Ra and 
Rb are the positions of their centers of mass and 



LA>niax{||r'-RA||} 

r' GA 



Lb ^ max{||r' — R_b||} 



(24a) 



(24b) 



Having calculated an overestimate of di^i^ (t) we can 
evaluate an overestimate of dAB that we call dmax- 



niax{d-n 



lAiE 



(25) 



Using Eq. (|25|l we can easily find an efficient strategy to 
bracket the solution. In fact the following algorithm can 
be used: 

1. Evaluate the distances between all sites that may 
interact {di^ig(t)}i^ig at time t (starting the first 
time from ti). 

2. Choose a time increment At as follows: 

r dABjt) 



At 



ifdAB(i)>e/; 
r^, otherwise. 



(26) 



where the two arbitrary parameters and e/ sat- 
isfy Ed < £/ < min{LA, Lb}- 

3. Evaluate the distances at time t + At. 
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4. If for at least one pair of patches (za, is) we find 
that the product di^i^ (t + At)di^ij^ (t) < we have 
bracketed a solution. We then find the collision 
times and the collision points solving Eq. H22I) for 
all pairs. Choose the smallest collision time and 
terminate. 

5. if pairs of patches are such that < \di^ig{t + 
At) I < ed and < \dij^ig{t)\ < e^, for each 
of these pairs evaluate the distance dij^i^it + 
At/2), perform a quadratic interpolation of these 
3 points {t,d^^^^{t)), {t + At/2,d,^,^{t + At/2), 
[t + At,dij^ig{t + At) and find if the resulting 
parabolas have zeros. If yes refine the smallest zero 
solving again Eq. ^221 for all these pairs. 

6. Increment time by At, i.e. 

t^t + At (27) 

7. Go to step l,ift< t2. 

If two particles undergo a "grazing" collision, i.e. a 
collision where the modulus of the distance stays smaller 
than ed during the collision, the collision could not be 
located by the previous algorithm due to failure of the 
quadratic interpolation. We have chosen to work with 
Cd ~ 10~^. For such a choice, we have not observed graz- 
ing collisions during the simulation runs (which would 
appear in the conservation of the energy). 

The basic algorithm can be improved with simple op- 
timizations. For example one can calculate d™"^ as fol- 
lows: 

dZZ - IIVabII + Wu^aWU, + W^bWU, (28) 

where 

L,^^\\y,J -B.a\\ (29a) 



^ + ApAB/A^(rA - xc) X n (31c) 



wb ^ ws - ApAB/^^(rB - xc) X n (31d) 

where n is a unit vector perpendicular to both surfaces 
at the contact point xc, 7^, Ib are the moments of in- 
ertia of the two colliding sticky particles, mA, tub their 
masses and the quantity ApAB depends on the type of 
the collision. If we define 

Vc = (vA+WAx(xc-rA)-VB-WBx(xc-rB))-n (32) 

If the collision occurring between particles is an hard-core 
collision, one has: 

ApAB = -2v, (33) 

if the collision occurred between two spherical patches 
already bonded (i.e. prior the collision the distance be- 
tween the two sites is < (5 one has: 

ApAB = l^'^^" ^ 2uo/Mred 

I —Vc + \/vc — 2uo/Mred othcrwisc 

(34) 

where 

^Ved = m^^+mg^+I^^\\{rA-xc)xh\\+Ig^\\{rB-xc)xn\ 

(35) 

Finally if the collision occurs between two patches that 
are not bonded (i.e. the distance between the two sites 
is > (5 prior the collision) we have: 

ApAB = -Vc + V^c - 2uo/Mred (36) 



L,, = |lri3'-RB|l (29b) 

and if c?ab (^) > the time increment can be evaluated 
in the following optmized way: 

At = min{d,,,Jt)/dZZ} (30) 

D. Collision of two particles 

At the collision time, one has to evaluate the new ve- 
locities of centers of mass and the new angular velocities. 
If xc is the contact point then the velocities after the 
collision can be evaluated as follows: 

va^ m^^ApABn (31a) 
vs ^ vs - rrig^ApABn (31b) 



VIII. APPENDIX: EVALUATING THE 
PRESSURE 

A. Evaluating the pressure in the ED code 

We define the quantity: 

AAa^pit) = V I Vc.p{t')dt' (37) 
Jo 

where Vai3 is the molecular pressure tensor, 

N N N 

ip - R]p) (38) 

i=l i—1 j>i 

The sums in the previous expression involve components 
(denoted by greek letters), Vi, Ri and Fij, which are the 
velocity, the position of center of mass of i-th. particle 
(mass Mi) and the total force acting between particle i 
and j respectively. 
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In the presence of impulsive forces, the stress tensor 
defined in Eq. (|38|l is not weU defined, while the integral 
in Eq. H37(l is well defined. Consider the time interval 
{t,t + At). During this interval the quantity AAapit) 
will vary due to the collisions occurring between particles. 
The variation SA{t) of AA{t), is: 

N 
i 

where 5t is the time elapsed from last collision occurred 
in the system and JP, is the variation of momentum of 
particle i after the collision, i.e. 



contribution, equal to nksT] a positive HS contribution, 
which requires the evaluation of the hard sphere radial 
distribution function g^^ (r) at distance cr and a nega- 
tive contribution arising from the SW interaction, which 
requires both the evaluation of the H — LP radial dis- 
tribution function gu-Lpi'"') distance 5 as well as the 
evaluation of < Rn-Lpir) >■ For a pair of H and LP 
sites whose distance is r, the quantity Rh-lp is defined 
as the projection of the vector joining the centers of the 
two particles associated to the two sites along the di- 
rection of the unitary vector joining the two sites. The 
ensemble average < ... > is performed over all pairs of H 
and LP sites at relative distance r nl. 



SPia = ApABfla 



(39) 



where ApAB is the quantity defined in Eq. 

From AAafj (t) and AA^p {t -f At) the average pressure 
over the interval At can be evaluated as follows: 



-E 



AA^^{t + At) ~ AA^^(t) 
At 



(40) 



B. Evaluating P in MC 



The resulting expression for P is 



P = nfcBT(l + -^- (41) 



9us{^) - 8^(1 - e-i/^) < Ru-Lp{^) > g^^.M 



) 



In the analysis of MC configurations, pressure has been 
calculated as sum of three contributions. A trivial kinetic 
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